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I. INTRODUCTION 

The last decade saw the successful construction and op- 
eration of the first generation of large-scale gravitational- 
wave interferometric detectors, at the observatories of 
LIGO IJ in the United States, and of Virgo [2] and 
GEO 600 [3 in Europe. Ongoing improvements to the 
sensitivities of these detectors (e.g. [4 ), combined with 
the construction of new large-scale gravitational-wave in- 
terferometers in Japan [5 and potentially elsewhere, are 
widely anticipated to result in the first ground-based de- 
tection of gravitational waves within the next decade. 
From these first detections will follow new tests of the 
fundamental physics of gravity, as well as new tools for 
exploring the astrophysics of compact objects. 

Gravitational-wave pulsars are one class of sources 
which may be detected by the next generation of 
gravitational-wave interferometers. They are rapidly- 
rotating neutron stars which emit long-lived, narrow- 
band, quasi-sinusoidal gravitational waves, and are of- 
ten referred to as continuous or periodic sources. The 
dominant gravitational-wave emission is expected to be 
due to non-axisymmetric distortions of the neutron star; 
other possible emission mechanisms are unstable oscil- 
lation modes such as r-modes, and free precession due 
to misaligned rotation and symmetry axes; see [5] for a 
review. An isolated neutron star may have acquired a 
non-axisymmetric deformation during birth in a super- 
nova, although it is uncertain for how long such a defor- 
mation might be retained. For neutron stars in binary 
systems, non-axisymmetry may result from e.g. differen- 
tial heating from accreted matter leading to differential 
density gradients [7 , or from the confinement of the ac- 
creted matter by the star's magnetic field in a magnetic 
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mountain [8l[9j. 

Signals from gravitational-wave pulsars are being ac- 
tively searched for in data from the LIGO and Virgo de- 
tectors. Although none have been found to date, many 
of these searches have reached sensitivities where (albeit 
under optimistic conditions) such signals could be de- 
tected. The most recent search for gravitational waves 
from known radio- and X-ray pulsars ^lOj constrained 
the gravitational-wave power radiated by the Crab pul- 
sar to be less than 2% of the total power available from 
the loss of rotational energy. A search for gravitational 
waves from the Vela pulsar [TT] set energy constraints of 
35-45%, depending on assumptions about Vela's orien- 
tation. Upper limits on the amplitude of gravitational 
waves from the neutron star in the supernova remnant 
Cassiopeia A, set in below the level expected (op- 

timistically) from the total conversion of rotational en- 
ergy into gravitational waves, assuming that the (as-yet 
unknown) rotation period of the neutron star is within 
the searched frequency band. Recent searches for undis- 
covered neutron stars which may be radiating gravita- 
tional waves [131 114j have set gravitational-wave ampli- 
tude limits comparable to upper limits hypothesized for a 
population of such stars [15]. Searches for gravitational- 
wave pulsars in binary systems have so far focused on the 
most promising known target, the low-mass X-ray binary 

scorpius x-1 [iniin]. 

To best assess the prospects of future searches for 
gravitational-wave pulsars, it is important to be able 
to accurately estimate the sensitivity such searches can 
achieve. (What is meant here by sensitivity is defined in 
In particular, designing searches for con- 
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tinuous gravitational waves which cover a wide range of 
possible signal parameters (e.g. searches for undiscovered 
neutron stars) commonly requires constructing a hierar- 
chical pipeline comprised of different data analysis tech- 
niques, each with different trade-offs, such as better sen- 
sitivity but increased computational cost, or vice versa. 
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An accurate estimate of the overall sensitivity of such a 
pipeline is important, therefore, for identifying the opti- 
mal combination of its elements. 

Obtaining the sensitivity of a search targeting a single 
source, such as a known pulsar, is relatively straightfor- 
ward (e.g. [B]). The calculation becomes more difficult, 
however, for searches over wide signal parameter spaces, 
e.g. searches for undiscovered neutron stars, or searches 
targeting known neutron stars with unknown rotation 
periods. This difficulty has resulted in two different 
approaches to sensitivity estimation. Wide-parameter- 
space searches of LIGO and Virgo data [l2HT4l [TB I [T8ti22] 
set upper limits on gravitational-wave amplitude (which 
in turn characterize the sensitivity of the search) by per- 
forming Monte Carlo simulations, where the search is re- 
performed on computer-generated data containing simu- 
lated signals (see Section II C I . While Monte Carlo sim- 



sis techniques used to search for them (Section II B I , and 



the method by which search sensitivities are estimated 



ulations are appropriate for accurately computing the 
sensitivity of searches of real gravitational-wave detector 
data, which often contain e.g. non-Gaussian instrumen- 
tal noise artifacts, they are usually too computationally 
expensive to be useful for theoretical studies of the sen- 
sitivities of different data analysis techniques. 

Instead, theoretical studies, e.g. p5H5T] . commonly 
make certain assumptions about the distribution of the 
gravitational-wave signals being searched for, in order 



to simplify the sensitivity calculation (see Section III C I . 
These assumptions, however, result in a measure of sen- 
sitivity that, as shown in Section HID is quantitatively 



different from that arrived at using Monte Carlo simula- 
tions. To date, there has been little published work in the 
gravitational-wave literature on the discrepancy between 
these two approaches. Furthermore, an accurate, com- 
putationally cheap (i.e. suitable for theoretical studies) 
estimator of sensitivity, as obtained by Monte Carlo sim- 
ulations, has yet to be proposed (although combinations 
of numerical and analytic sensitivity estimation methods 
have been developed; see [32l[33]). It is these two issues 
that the present work intends to address. 

Section [ll] of this paper presents an overview of 
gravitational-wave pulsar searches, and describes the 
most common method by which their sensitivities are 
estimated. Section IIIII contains the main result of this 
paper: an analytic expression which may be used to 
quickly and accurately estimate the sensitivity of wide- 
parameter-space searches for gravitational-wave pulsars. 



Section IV verifies the accuracy of the analytic sensitiv- 



ity estimator, and Section [V] assesses the validity of some 
assumptions that were made during its derivation. Sec- 
tion |Vl] discusses the results presented in this paper, as 
well as possible avenues for future research. 



II. GRAVITATIONAL- WAVE PULSAR 
SEARCHES 

This section is an overview of the signal model of 
gravitational- wave pulsars (Section |lI A ), the data analy- 



(Section II C). See also [B] for a review of gravitational- 
wave pulsar data analysis, and [5T for an overview of the 
current data-analysis activities of the LIGO and Virgo 
scientific collaborations. 



A. Signal model 

The signal from a gravitational-wave pulsar is written 
as a time series h{t) of the dimensionless strain ampli- 
tude h, which for a ground-based interferometric detec- 
tor is proportional to the differential change in the length 
of its arms (which are assumed to be much shorter than 
the gravitational wave-length). We assume that the sig- 
nal contains only a single frequency component, although 
it is possible for it to contain multiple frequency compo- 
nents arising from free precession |35j . Following ^23^ , the 
time series h{t) may be written as the summed products 
of four amplitudes Ai and four time-dependent functions 

h^it)■. 



h{t) = ^AMt) 



(2.1) 



The amplitudes Ai are related to the four amplitude pa- 
rameters of the signal; its overall strain amplitude /iq; 
its initial phase (/)o; the inclination angle t between the 
neutron star angular momentum and wave propagation 
vectors; and the polarization angle ip, which fixes the ori- 
entation of the neutron star about the wave propagation 
vector. The functions hi{t) are functions of the signal's 
remaining phase parameters: its sky position, given by 
its right ascension a and declination (5; and its frequency 
evolution, given by an initial frequency /, and frequency 
time-derivatives or spindowns /, /, f^^\ etc. The num- 
ber of spindowns required generally depends on the age 
of the sources being targeted [5BH55] . 

The time series h{t) may also be written in a form 
which illustrates the two polarizations, plus and cross, of 
a gravitational wave: 

h{t) = A+F+(i)cos$(t) +^x-Fx(t)sin$(t) , (2.2) 

where A^ and A-^ are the amplitudes of their respective 
polarizations, and ^(t) is the signal phase. The antenna- 
pattern functions F-^-{t) and Fx(t) give the response of 
the detector to each polarization, and are modulated by 
the sidereal motion of the Earth. Expressions for (t) 
and Fx (t) are given in [531 , and in Appendix |Aj 

The signal-to-noise ratio (SNR) of a signal, p, is found 
by integrating h(t) over the observation time T, which 
gives: 



nT/2 



a. 



(2.3) 



where (F^)^ = {1/T) J^^^^^dt F{t) 
of F+(t)2 (similarly for {F^ 



^ is the time average 
- j4+//io (similarly 
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for Ox), and Sh is the one-sided detector noise power 
spectral density. The SNR is independent of the signal's 
phase modulation; when T is longer than several days, 
the SNR's dependence on the signal's amplitude modula- 
tion also vanishes, and becomes a linear function of T. 
Expressions for {F^)t and {F^)t are given in Appendix [a] 
see also |23| for expressions for p^. 



B. Search techniques 

Gravitational-wave pulsars are very weak sources; the 
amplitude of their signals is likely to be several orders of 
magnitude smaller than the noise amplitude of current- 
and even next-generation interferometric detectors. Nev- 
ertheless, gravitational-wave pulsar signals can be recov- 
ered using the well-known technique of matched filtering, 
where the data are correlated against a template which 
models the signals' amplitude and phase evolution over 
time. Matched filtering was first applied to the detec- 
tion of gravitational- wave pulsars in |23| . and extended to 
multiple detectors in [13] . The detection statistic derived 
in these papers, known as the J^-statistic, maximizes the 
signal SNR over the four amplitude parameters, but re- 
quires values to be chosen for the phase parameters. 

A search for gravitational-wave pulsars therefore con- 
sists of performing matched filtering against a bank of 
templates, whose phase parameters are chosen to cover 
the parameter space of interest, e.g. over the whole sky 
and a broad range of frequencies for a search for unknown 
neutron stars. It is almost certain, however, that any sig- 
nal in the data will possess parameters which are different 
from any one of the searched templates; consequentially, 
no template will perfectly match the signal, and the sig- 
nal SNR will be degraded. While some loss in SNR is 
unavoidable, template banks are constructed such that 
the fractional loss in SNR, also known as the mismatch, 
can never be greater than some prescribed maximum. To 
accomplish this, a metric is often used to determine how 
closely the templates must be spaced in each parameter 
[361 144j . How to construct a bank which minimizes the 
number of templates is known in theory [45^ , but is often 
difficult to accomplish in practice. 

Unfortunately, the number of templates which must 
be matched filtered increases rapidly with the length of 
data being analyzed. If T denotes the time-span of the 
analyzed data, the number of templates which must be 
placed in each parameter dimension scales as follows: 
for the sky position (a, J), T for frequency, and 
rpk+i each spindown parameter Z'*^^. For example, 
the number of templates for an all-sky search requiring 
one spindown scales as T'^^. In contrast, the sensitiv- 
ity achievable by matched filtering increases only as T^/^ 
(assuming that the data are contiguous in time). Due to 
the rapid increase in the number of required templates, 
matched filtering quickly becomes too computationally 
expensive for searching long data sets and large parame- 
ter spaces [23 ■ 



The solution is to resort to a hierarchical pipeline, 
where typically the data are broken into short segments, 
each of which are matched filtered separately. The re- 
sults from each segment are then combined using semi- 
coherent analysis methods, which resemble matched fil- 
tering but do not require full amplitude and phase con- 
sistency of the signal template between data segments: 
typically, only the derivative of the phase (i.e. the fre- 
quency) must be consistent. (To distinguish it from semi- 
coherent methods, matched filtering is also referred to as 
coherent matched filtering.) For fixed T, and assuming 
no limits on computational cost, semi-coherent methods 
are less sensitive than coherent matched filtering: their 
sensitivity scales roug hly as Ns^^Ts^'^, where and 
are the number and time-span of each segment, while 
a coherent search (assuming contiguous data) scales as 
~ Ns^^Ts^^. On the other hand, wide-parameter- 
space gravitational-wave pulsar searches are almost al- 
ways computationally limited, due to the large param- 
eter spaces which must be searched. The number of 
templates, and hence the computational cost, of semi- 
coherent methods scale with Tg <C T, instead of T, mak- 
ing them computationally cheaper than a fully-coherent 
search. This in turn permits a semi-coherent search to 
use more data, improving its sensitivity, while remaining 
computationally tractable. Some hierarchical searches, 
e.g. [201 [22], use Tg ~ 30-minute segments, so that 
the coherent matched-filtering step closely resembles the 
computation of a power spectrum; other searches, e.g. 
[131 m] use segments of Tg ~ 1 day. 

Examples of semi-coherent methods are the Stack- 
Slide [MlUZlllH], Hough f^^,"^, PowerFlux [M], cross- 
correlation [37], and global correlation transform [iS] 
methods. Recently, methods which blend together as- 
pects of semi- and fully-coherent methods have been de- 
veloped ^49. - i51j . The problem of how to optimise hierar- 
chical searches is studied in [371 EI] ■ 

C. Search sensitivity 

The sensitivity of a wide-parameter-space search for 
gravitational- wave pulsars has traditionally been charac- 
terized by the method presented in this section; two al- 
ternative methods of estimating sensitivity are discussed 
in Section [VH 

The method follows from the canonical framework for 
statistical hypothesis testing formulated in [321 US! ■ It is 
commonly referred to as the frequentist method in the 
gravitational-wave literature. In essence, the method 
provides the answer to the following question: if there 
were a population of gravitational-wave pulsar signals 
present in the searched data, each with the same ampli- 
tude ho, how large would ho need to be before we would 
be confident of detecting a very large fraction of them, 
e.g. 95%? 

In order to make the above question more precise, we 
must first define what is meant here by detection. This 
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is complicated by the unavoidable fact that the output 
of any real gravitational-wave detector contains noise, in 
addition to any signal. When noisy data is analyzed, 
e.g. using matched filtering, the results of the analysis 
may be subject to false alarms and false dismissals. A 
false alarm is when the results of the analysis falsely in- 
dicate the presence of a signal, which is instead simply 
due to spurious noise fiuctuations; a false dismissal is 
when a real signal present in the data is sufficiently cor- 
rupted by noise that it fails to be identified as a signal 
in the analysis results. We define a detection statistic to 
be a number quantifying a single search result, e.g. a 
correlation of the data against a single signal template. 
Typically, the detection statistic increases with the prob- 
ability that the data contains a signal which matches the 
template. As an example, the J^-statistic represents the 
result of a matched-filtering analysis maximized over sig- 
nal amplitude parameters. We say that a signal has been 
detected^ if the value of its detection statistic s exceeds 
some threshold S[^ 

Suppose that we perform a wide-parameter-space 
search of data which is known to contain no signal, e.g. 
the data may comprise only computer-generated Gaus- 
sian noise. Because the data used to compute the detec- 
tion statistic are combined differently for each template, 
the values of s returned by the search will not be identi- 
cal, but will instead follow a certain probability distribu- 
tion. The probability of a detection in this data, i.e. the 
probability of a false alarm, denoted pfa, is the proba- 
bility of the detection statistic s exceeding the threshold 
Sfa, under the assumption that no signal is present in the 
data: 



Pfa 



p{s > Sfa I no signal in data) 



(2.4) 



where p{A\B) denotes the probability of the statement A 
being true, given that we already know that B is true. 

We now perform the same search using data which is 
known to contain a single signal, e.g. a simulated signal 
added to computer-generated Gaussian noise. Assuming 
that the majority of the searched templates are insen- 
sitive to the signal (which would be true of a properly 
constructed template bank), the majority of the returned 
values of s will follow the same probability distribution 
as before. For templates whose parameters are close to 
those of the signal, however, the values of s will follow a 
different probability distribution, which instead assumes 
that a signal is present. Assume that the SNR of the sig- 



^ It is important to note that, when performing a search of real 
gravitational-wave detector data, any signal thus detected would 
never be automatically claimed as a genuine gravitational-wave 
signal; extensive follow-up investigations would first be per- 
formed to e.g. consider possible contamination by instrumental 
artifacts. 

^ The hierarchical searches in 1131 1211 employ a slightly differ- 
ent definition of detection: the detection statistic must exceed 
a threshold in a given number of data segments. 



nal is large enough that these templates, close to the sig- 
nal, are distinguishable from the remaining search tem- 
plates; this assumption is satisfied by requiring a low false 
dismissal probability. The probability of this particular 
signal not being detected, i.e. the probability of its false 
dismissal, denoted Pfd(p), is the probability of s, in the 
neighborhood of the signal, falling below the threshold 
Sfa, under the assumption that a signal with parameters 
p is present in the data: 

Pid{p) = p{s < Sfa I signal(param. p) in data) . (2.5) 

Note that we assume here that the search template is 
perfectly matched to the signal, whereas in reality there 
will be some mismatch between them. (The loss of sen- 
sitivity due to mismatch between template and signal 
is considered in Section |V C ) We repeat the search for 
a large number of signals, each with different p drawn 
from a chosen distribution of signals parameters, p{p). 
The overall false dismissal probability, denoted pi^ , is the 
average false dismissal probability obtained from each of 
the sampled signals: 



Pfd = (pfd(p))^ = J ■ 



dpp{p)pM{p). (2.6) 



Equation (2.6) is computed using Monte Carlo integra- 
tion, i.e. by the repeated computation of Eq. (2.51, with 



different parameters p drawn from the distribution p{p) . 

We now return to the question stated at the begin- 
ning of this section: how large would the amplitude ho 
of a population of gravitational wave signals need to be 
in order fo r th em to be detected, e.g. 95% of the time? 
First, Eq. (2.4) is solved for the threshold Sfa that would 
result in a desired false alarm probability pfa. Then, 



given a target false dismissal probability ptd, Eq. (2.6) is 
solved for the overall strain amplitude Hq. Suppose h^^'^" 
is the value of ho which solves Eq. (2.6 1 for ptd = 5%; 



then, if a population of signals have constant amplitudes 
Hq^'^", then a fraction 1 — pfd = 95% of them will be de- 
tected. The amplitude hg^'^" thus characterizes the am- 
plitude of signals a particular search method is able to 
confidently detect, and hence gives a useful measure of 
the search's sensitivity. Where a search finds no credible 
gravitational- wave signal, this sensitivity is interpreted as 
an upper limit, with confidence 1 —ptd, on the amplitude 
of signals present in the searched data. 

The above procedure is used to set upper limits on 
the amplitude of gravitational waves for wide-parameter- 
space searches of LIGO and Virgo data (see Section |l] 
for references). The search parameter space is typically 
partitioned into small frequency bands, and upper limits 
are set separately for each frequency band. One impor- 
tant difference to the procedure described in this section 
is that Sfa is determined by the largest value of s re- 
turned by the search (after instrumental noise artifacts 
have been removed); an effective false alarm probability 
can then be determined from Sfa using Eq. (2.4). 
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III. ANALYTIC SENSITIVITY ESTIMATION 

Theoretical studies of gravitational-wave pulsar search 
pipelines require a method of efficiently and accurately 
estimating the sensitivity achievable by such searches. 



given by Eq. (2.3) for some appropriate value of Sh- For 



The procedure described in Section If C is generally un- 
suited to this task, due to the computational expense 
of repeatedly generating and searching data containing 
simulated signals. In this section, we construct an ana- 
lytic expression which accurately reproduces sensitivity 
estimates computed using the frequentist method, using 
the procedure described in Section |II C[ but is simpler 
to implement and computationally cheaper to calculate. 
Expressions for the false alarm and false dismissal prob- 
abilities are presented in Sections |III A| and |III B| re- 
spectively. Section |III C| presents a commonly used, but 
inaccurate, analytic estimator of the sensitivity, and Sec- 
tion [1110] presents a new, more accurate expression. 

We restrict our attention to detection statistics s which 
follow distributions (e.g. [Mj)- This imphes that s 
can be written as the sum of squares of some number 
of normally-distributed random variables. Examples of 
such statistics are the StackSlide power j28], and the J"- 
statistic [231 [57]. Our prototypical search is a single-stage 
hierarchical search of Ng data segments, each of which 
span a time Tg. (A fully-coherent search is then given 
by the special case Ng = 1.) A coherent analysis is per- 
formed for each segment i, returning detection statistics 
Si , which are then summed using a semi-coherent method 
to attain the final detection statistic s — X^i^i ^i- 

In the absence of a signal, s is distributed according to 
a central distribution, which takes a single parameter: 
the number of degrees of freedom of the statistic. Since 
s is the sum of the Ng values s;, its number of degrees 
of freedom is given by NgV, where each of the Si has v 
degrees of freedom. A detection statistic derived from the 
power of a signal (e.g. the StackShde power) has = 2; 
the J^-statistic has v = A. We denote that s is distributed 
according to a central distribution with NgV degrees 
of freedom by s '-^ x^{Ngi', 0). 

In the presence of a signal, s is distributed according to 
a non-central distribution, which take two parameters: 
the number of degrees of freedom NgV as before, and the 
non-centrality parameter 



Ngp' = hlT{al{Fl)t^ 



JVdet. 

^ ^(5,)-.i. (3.1) 



The non-ccntrality parameter is the accumulated SNR of 
a signal analyzed over Ng data segments from A^dct. de- 
tectors, assuming perfect match between signal and tem- 
plate. The noise power spectral density of the ith data 
segment from the jth detector is labeled {Shjij- It is con- 
venient to define NgS^^ = T.f^^T.'^lT iSh)~\ where- 
upon we recover the expression for given in Eq. (2.3). 



For simplicity, therefore, we can ignore the summation 
of {Sh)ij over segments and detectors, and take to be 



a real detector, Sh also varies as a function of frequency; 
we therefore assume that we are considering the sensitiv- 
ity of a search over a frequency band sufficiently narrow, 
so that Sh can be assumed constant. 

We denote that s is distributed according to a non- 
central distribution with NgV degrees of freedom and 
non-centrality parameter Ngp'^ by s ~ x^{NgV,Ngp'^). 
Note that, as in Section jll C| we assume that the search 



template is perfectly matched to the signal, and delay 
considering the sensitivity lost due to template-signal 
mismatch until Section [V CI 



A. False alarm probability 

We start by determining the threshold Sfa on the detec- 
tion statistic required to give a certain false alarm prob- 
ability Pfa. 

Assuming that no signal is present, the probability that 
a single value of the detection statistic s falls below sta 
is p(s < Sfa|s ~ x^lNgV.Q^. Suppose that the search 
returns values of s, i.e. A^t template^ are searched. 
We assume that each value of s is statistically indepen- 
dent, i.e. that the joint probability of obtaining any two 
values si and S2 is just the product of the probabilities 
of obtaining si and S2 individually. (The validity of this 
assumption is examined in Section VA ) The probabil- 
ity that A't values of s fall below Sfa is then given by A't 
multiples of the single- value probability: p{---)^*- Its 
complement, 1 — p(---)^S the probability that in a 
search of A't templates, one or more values of s will be 
returned above the threshold, i.e. the probability of a 
false alarm, pfa. In short, we have 



Pfa = 1 - [p{s < Sfa|s ~ x'lNsiy.O])]^' . (3.2) 

Simple re-arrangement gives 

(1 -Pfa)^/^' = p{s < sfa|s ~ x^[Ngiy,0]) . (3.3) 

Since typically either A't = 1 (e.g. for a search f or a 
known pulsar) or A't ^ 1, the left-hand side of Eq. (3.3) 



can be replaced with the first-order binomial expansion 
1 — Pfa/Nt- We now have 



PfaM = P{s > Sfa|s ~ X^[Ngiy, 0]) 



(3.4) 



Note that Pfa./Nt can be interpreted as the false alarm 
probability for a single template out of the A't templates 
searched. 

We must now solve Eq. ( |3.4| for Sfa, i.e. we must com- 
pute the inverse function of the central distribution. 



^ Note that we are referring to the number of templates searched 
in the semi-coherent stage of a hierarchical pipeline, which is 
distinct from the number of templates searched in each individual 
data segment; see e.g. |27| . 
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FIG. 1. Normalized false alarm threshold function 
of Pfa/A^t and A'a, with v = 4 (thick lines). Its limiting values 
for large A'a, given by the normal distribution approximation, 
are plotted as thin vertical lines for Zfa = 3-7. 



For small values of Pfa/Nt (typically, pta is chosen to be 
1%, and Nt ^ 1), an analytic expression for Sfa derived 
from 55 is given in Appendix [Bj We define the normal- 
ized false alarm threshold Zfa to be the difference between 
Sfa and the mean of a central distribution with NgV 
degrees of freedom, in units of the distribution's standard 
deviation: 



Sfa - NsV 



(3.5) 



In the limit of large TVs, -Zfa ~ ■\/2erfc^^(2pfa/iVt) and is 
independent of A's (see Appendix [b|). Note that this lim- 
iting value for Zfa is identical to the expression obtained if 
the central distribution is approximated by a normal 
distribution, as is commonly done (e.g. in [H]). Figure [l] 
plots Zfa as a function of Pfn/Nt and Ng, and illustrates 
its convergence to the normal distribution approximation 
for large A's. 



B. False dismissal probability 

Having determined the threshold Sfa appropriate for a 
desired false alarm probability pfa, we now attempt to 
solve Eq. (2.6) for some quantity which usefully charac- 

The quantity most 



terizes the sensitivity of the search, 
often used for this pu rpose has been the dimensionless 
expression ho y/Tg/Sh, which is sometimes referred to as 
the statistical factor [301 138j . The statistical factor quan- 
tifies the sensitivity of a search to a population of signals 
of amplitude ho, relative to the performance of the de- 
tector given by its noise power spectral density Sh, and 
the length of (coherently) analyzed data T^. 

In this paper, we propose instead to use the root-mean- 
square SNR, \/{f?), to quantify sensitivity. H ere, () 
denotes averaging over a, sin 5, ip [see Eq. (A9)], and 



-0 
Si, 
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FIG. 2. False dismissal probability pfd as a function of 
\/(^, given by Eq. (|33|, with u ^ 4, and for pfa/iVt = 0.01 
(sohd) and lO'^^ (dashed), and iV^ = 1 (left) and 10"* (right). 



of averaging implies that the population of gravitational- 
wave pulsars being searched for are isotropically dis- 
tributed over the sky, and are isotropically oriented. Us- 
ing Eq. (2.3), we have 



a, sin 5,ip.t 



+ {a\)^{Fl)a.s\nS,-iP,t) ■ (3.6) 



The averages of F'^ and over a, sin (5, and ij: are given 
by Eq. ( All[ ) (where we assume C, — 7r/2); after this av- 
eraging Fj^ and are independent of time. We assume 
a signal generated by a non-axisymmetrically deformed 
neutron star, for which (e.g. [23l 135]): 



i+e 



(«+>« 



(4>« 
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Finally we have 



(3.7a) 
(3.7b) 

(3.8) 



^ = cosi, given by {p^)^ = (1/2) /j'^^df xhis choice 



i.e. {p^) is directly proportional to the statistical fac- 
tor. Unlike the statistical factor, however, \f{fP^ relates 
directly to a property of the population of signals being 
searched for (i.e. their mean SNR), and hence is a more 
directly physical quantity. It also has a clearer interpreta- 
tion as a measure of sensitivity: for example, to improve 
the sensitivity of a search, we must make the search able 
to detect signals with weaker SNR (at the same false 
alarm and dismissal probabilities), and hence we must 
lower the mean SNR, i.e. \/~{p^, of the population of 
signals which the search can detect. It is convenient to 
write p in terms of and a factor i?, defined such 

that p — \J {p'^)R, which implies {R^) — 1. 

Assuming that a signal with parameters p is present, 
the probability that the detection statistic s (in the neigh- 
borhood of the signal, as discussed in Section II C ) falls 
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below Sfa is Pfd(p) = p{s < Sfa|s x^Wsi',Ns{p'^}R'^]). 
Since ho, Sh, and Tg are taken to be constants, (p^) is 
also a constant. The only quantity which depends on the 
signals parameters is therefore R, which is a function of 
p — {a,S,tp,£^). By Eq. (2.6|, the overall false dismissal 
probability pf^ is 



p(s < Sfa S ' 



a, sin S,'ip,$^ 



(3.9) 



3.9 1 is plotted in Fig. [2] for different choices of 



Equation 
Pfa,/Nt and „ 

To proceed, we must now solve Eq. (3.9) for (p^). Un- 
fortunately, an analytic solution is difficult to obtain, due 
to the complicated dependence of R on the signal param- 
eters. Two approaches to solving Eq. (3.91 analytically 
are presented in Sections |III C| and |III D 
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FIG. 3. Contours of contributions of the terms T (see the 
text) to the second factor of Eq. (3.131: \/Ti/ 'Y^^ Ti (solid 
contours), ^JTij^^Ti (dashed contours), and 's/Ts/Yji'Ti 
(dotted contours). 



C. Sensitivity to constant-SNR signal populations 



To solve Eq. (3.91, it is common to assume that ev- 
ery signal, in the population of signals being searched 
for, has the same SNR, which we denote by p. An al- 
ternative interpretation of this approximation is that the 
population of signals can be replaced by a single signal, 
whose SNR is p = {p^)- While this assumption does 
not accurately model a physically reasonable population 



of gravitational- wave pulsar signals (see Section 111 D ) , it 
does allow Eq. (3.9 1 to be readily solved for p: because 



the signals all have the same SNR, no averaging over sig- 
nal parameters is required, and {p^)R^ is simply replaced 
by . The accuracy of this approximation is examined 
in Section Imp] 

Another common simplification is to approximate the 
non-central distribution by a normal distribution 
Af{fJ., (j) with the same mean and standard deviation, 
which are /i = Ns{i' + p^) and a = ^/2N^^Jj7+2]P) re- 
spectively. This approximation introduces an error in p 
of < 2.5% (at pfa/iVt = 0.01, Ns = 1) which decreases 
with decreasing Pfa/N^ and increasing Ng. With this ap- 
proximation, Eq. (|3.9|) reduces to 



(3.10) 
(3.11) 



Pfd wp(s < Sfa|s 7V[/i,a-]) 

^ 1 ^^^^ / N,p^ - zi^ym^ ' 

2 % 2ViV.(i. + 2p2) ^ 



where we have substituted Sfa with the normalized 



threshold Zfa defined by Eq. (3.5), and crfc is the com- 



plementary error function. The solution to this equation 
is 



,2 r 



l + (l + Q)- 



(3.12) 



where Q_= {N^u + Zf^V8N,i^)/{2q^), q = v^erfc-^ 2pM, 
and erfc~i is the inverse complementary error function. 



For pfa/A^t < 0.01, pfd > 0.05, > 1, a nd > 2, Q > 5, 
and it is reasonable to approximate ^1 + Q with y/Q. 
This permits further simplification to 



2i/ 




(3.13) 



- —1/4 

The first factor of p scales with Ns , which is a well- 
known property of hierarchical searches [551 [25] • The 
second factor contains, inside the square root, a con- 
stant term, Ti = Zfa, a term which scales approximately 
with N^^^'^, Ti ^ q{l + Zf^y/S/N~uYl'^ , and a term 
which scales with Ns , Ts, ~ q^ \/2/NsV. These addi- 
tional terms appear because we do not employ the weak- 
signal approximation of [26j . which Taylor-expands pfd 
[Eq. (3.11)] to first order in p. This approximation is not 



valid here because we are interested in small false dis- 
missal probabilities, which implies that a large fraction 
of signals will be strong, i.e. enough to cross the detec- 
tion threshold. Each of the terms T contribute to the 
value of p over the ranges pia/Nt < 0.01 and iVg > 1, as 
shown in Fig.jsj over the plotted ranges ofpfa/Nt and A''s, 
the contributions of each term to the sum is 70- 

90% for /7T, 40-60% for V%, and < 40% for ^/%. For 
small pta/A*t and large Ng, \/7T dominates the sum, and p 



begins to scale purely with N. 



-1/4 



A general power-law 



scaling of p with Ns is utilised in [31 . 



D. Sensitivity to isotropically-distributed signals 



This section presents a new, more accurate sensitivity 



estimator than that presented in Section 111 C The im- 
portant difference is that we will no longer assume that 
all signals being searched for have the same SNR, and 
instead perform the correct averaging of ptd (p) as given 
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FIG. 4. The distribution of -Ri^o, calculated analytically 
using Eq. (3.151 (thin black line), and numerically (thick 
gray line). The black dashed line plots the distribution of 



so,mism 



^iso(l ~ m)i where the mismatch ^ is drawn from 
a distribution appropriate for a lattice-based template bank 
(see Section V C I . 
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FIG. 5. False dismissal probability pfd(C) oi a single signal 
observed in an isotropically-sensitive detector network, as a 
function of the cosine of the signal's inclination angle, ^ = 
cost. The exact result of Eq. (3.171 (solid line) is compared 
against the Gaussian function approximation of Eq. ( |3.18[ ) 
(dashed hne), for pfd = 1%, v = A, and: pfa/iVt = 0.01, iVs = 
1 (left), and piJN, = IQ'^^ , N, = 10* (right). 



in Eq. (3.9). This implies that we are searching for a 
population of signals isotropically distributed in sky po- 
sition (a, i5) and orientation {tp,£,) parameters. Unlike 
the assumption of constant SNR, this is a more physi- 
cally reasonable assumption. The observed distribution 
of nearby (millisecond) radio pulsars is roughly isotropic 
in the sky [56], and the angular momentum of neutron 
stars (which determines ip and ^) is not expected to have 
a preferred direction. 

We also assume that we are searching data from a net- 
work comprising a large number of gravitational-wave de- 
tectors, evenly distributed over the Earth, such the net- 
work is isotropically sensitive to gravitational wave arriv- 
ing from all directions. (The validity of this assumption 
is investigated in Section [VB ) We can therefore approx- 



imate by [see Eqs. A9 and AlO 



-^iso(0 — )'I>a,sinA,7 — {R )a,sinS,^ 



16 



(3.14) 



Substituting into Eq. (3.9 1, we see that we now need only 
to average Pfd(p) over ^. Figure |4] plots the probability 
distribution of i?? ^ assuming a uniform distribution in ^, 
given analytically by 



ISO/ 
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(50 + 20i?2 J5 - 7(50 + 20i?L) , (3-15) 



where Rf^^ ranges from 5/16 (^ = 0) to 5/2 (|^| = 1). 
Note that the most probable values of i?? ^ are those from 
linearly-polarized signals (^ = 0), with a rapid fall-off to- 
wards circularly-polarized (|^| = 1) signals. It is clear 
from this plot that the assumption that all signals have 
the same SNR (i.e. that i??^ = 1), as assumed in Sec- 
tion |lirc] is not a reasonable one. 

We next approximate the non-central distribution 



by a normal distribution Af, as per Section III C Here, 



the approximation introduces a smaller error of < 0.5% 
(at Pfn/Nt = 0.01, A's = 1) which also decreases with 
decreasing pfa/A't and increasing Ng. Equation (3.9) now 
reads 



Pid 



where 



Pfd(C) 



■ erfc 



(3.16) 



(3.17) 



and we define p = \J (p^). 

We find that pfd (C) is well-approximated by a Gaussian 
function: 



(3.18) 



for suitable choices of A and B. This is illustrated in 
Fig. [5) where we plot Eq. ( |3.17[ ) against Eq. ( |3.18| for 

two different choices of Pfn/Nt and Ng. We note that it is 



important only that the integral of Eq. (3.18) accurately 
approximate that of Eq. (|3.17 ); thus, the slight underesti- 
mation of j5td(^) by Eq. (3.18) at ^ w is partly offset by 
its overestimation at |f | > 0.5. We note that Pfd(C) ~ 
at 1^1 = 1; numerical investigations confirm that this 
property holds true for Pfa/A^t 0.01 and Ng > 1. 



We can therefore replace the integral of Eq. (3.18) over 
— 1 < ^ < 1 with one over —00 < ^ < 00, since the in- 
tegral over 1^1 > 1 contri butes little to the value of pfd. 
The integral of Eq. (3.18) then simplifies to 




(3.19) 



We choose A by setting ^ = in Eq. (3.18), and obtain 



A = 



2^N,{v + 2pRl) ^ 



rfc 



(3.20) 
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where i?Q = Rfs^iC = 0)- Next, we choose > to be 
the value of ^ such that Pfd(C = £,i) [Eq. (3.17)] equals 



Note the similarity of the right-hand sides of Eqs. (3.26) 
and (3.28). We now define the effective false dismissal 



the targe t fals e dis missal probability pfd. By equat- 
ing Eqs. (3.11) and (3.17), we deduce that satisfies 



p'^Rfsoi^ = Ci) = P^ and is given by 



probability pj^, normalized false alarm threshold z^^, and 
number of segments N!.: 



6 = 




- 3 



Pfd = Ptd 



(2pfdS) 



Zf^r, Ni^N,r\ (3.29) 



(3.21) 



Using these quantities, we equate Eqs. (3.26) and (3.28), 
obtaining 



Substituting into Eq. ( |3.18[ ) gives the solution for 

B = ^ In ( — 1 . (3.22) 



^1 \PtdJ 



erfc^^2pj, 



'fd ■ 



(3.30) 



Finally, by noting that Eq. (3.30) is similar in form to 



Combining Eqs. (pgl), (pOl), and and taking Eq. (|3Jlj), it follows that its solution for p is given by 



logarithms, we have 
In ^ ^P'^ 



^ , =\nA In 



In 



A 
Pfd 



(3.23) 



The terms on the right-hand side of Eq. ( 3.23 ) may be 
Taylor-expanded to second order in x = erfc^2A: 



lnA = -ln2- -^x - ^x^ + 0{x^), 
a/tt tt 



In 



In 



A 

Pid 



= hi VI ln2p; 



1 



'fd| 



/tt In 2p 



•fd 



'-n^2-' + 0{x^)- 
7r(ln2pfd)2 



(3.24) 



(3.25) 



where the second expansion is valid only for pfd < 1/2. 
For pfa/iVt < 0.01, Pfd > 0.05, > 1, and v > 2, these 
expansions introduce errors of < 1% and < 16% respec- 
tively. Substituting Eqs. (|3.24[) and (|3.25| into Eq. (|3.23[), 



and using Eq. (3.20), we solve for 

^ N.p^Rl - zi,V2K^ 
2^N,{u + 2p-^Rl) 

= ^r-i (Vl-2(l-A)ln(2pfdS)-l 



(3.26) 



where 



F 
A 



1 



In 2pfd 
1 



1 + 21n2pfd 
2 



+ 



l + 21n2pfd (l + 21n2pfd)2 
2 



In 2p, 



'fd 



(3.27a) 
(3.27b) 

(3.27c) 



Now, suppose that p'^^ — (1/2) erfcy, and Taylor-expand 
Inpf^ to second order in y, as in Eq. (3.24); solving for y 



gives an expression (valid only for p[^ < \/e/2) for 

y = erfc-i 2pJd - ^ Ul-2\n2p[^ - l) . (3.28) 



analogy to Eq. (3.131 



1 

Rq 



2v 




(3.31) 



where q' = v^erfc ^2pL, 



Evaluation of Eq. (3.31 ) is complicated by the fact that 



its right-hand side is itself a function of p, through p[^ 
[Eq. ( [3^ ], S [Eg . (|3]27c])], and [Eq. pT2l] )]. Never- 
theless, Eq. ( |3.31 1 may be iteratively solved for p using 
the following scheme. First, a reasonable initial guess. 
Pa is chosen: an appropriate choice is /5o ~ 1.4p (see 
Fig. [7|. Next, po is used to compute an updated value, 
pi, by substituting into Eq. (3.31): pi = p{po)- There- 



after, new values of p are obtained using the mean of the 
previous two values, i.e. the nth value of p is 



Pn P 



Pn-l + Pn-2 



(3.32) 



Using this scheme, the sequence of values {pn} reliably 
converges to an accurate value of p: typically, 20-80 it- 
erations are required to achieve sufficient accuracy. The 



use of the mean of the previous two p in Eq. (3.32) sup- 
presses divergent oscillations in the sequence {/5„}. 

Figure |6] plots p as a function of Pfa/Nt and TVs, with 
Pfd = 0.1, for a detection statistic with v = i degrees 
of freedom (e.g. the J^-statistic) . For a single-template 
search at 1% false alarm and 10% false dismissal proba- 
bilities, we see that pfa = 0.01, Nt = Ng = 1, and p — 6.3, 
which should be interpreted as the average sensitivity of 
a collection of single-template searches which cover the 
parameter space uniformly in sky position (a, 6) and ori- 
entation (ipT^). This value of p is equivalent to a statis- 
tical factor of 15.7. This differs from the often-quoted 
(e.g. [IH]) statistical factor for a single-template search 
of 11.4, because that statistical factor is calculated as- 
suming a signal with an average SNR, i.e. it is calculated 
from p instead of p. For the search for Cassiopeia A pre- 
sented in HH, pfa = 0.01, Pfd = 0.05, Nt - 1.8 x 10^° [35], 
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FIG. 6. Contours of p as a function of pfa/A^t and A^s, with pfd = 0.1 and v — 4. Thick sohd contours are in units of 1.0, thin 
solid contours are in units of 0.5, and thin dashed contours are in units of 0.25. 
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FIG. 7. Contours of the ratio of p/p as a function of pta/Nt 
and A'a, with pfd = 0.1 and v = A. Unlabeled thin contours 
are in units of 0.01. 



Ns — 1, and p « 14.5, which implies a statistical factor 
of ^ 36 consistent with that quoted in [33l |38] . 

Figure [7] plots the ratio of p to p, as a function of 
Pfa./Nt and Ng. Relative to p, we see that p underesti- 
mates the mean SNR detectable by a search by ^ 29±5%, 
i.e. it overestimates the search's sensitivity by the same 
amount. For comparison, the typically amplitude cali- 



bration error of the LIGO detectors is ^ 10% [57], and 
the ratio of the best upper limits on gravitational waves 
from Cassiopeia A to the indirect limits from energy con- 
servation is ^ 60% [12J. Thus, an error of ~ 30% in 
estimating a search's sensitivity is a significant discrep- 
ancy. Note, however, that the change in p/p is small over 
the ranges of pta/Nt and Ng plotted in Fig. [7j We con- 
clude from this that, while p does not predict the correct 
sensitivity, it does capture the correct scaling of sensitiv- 
ity with the false alarm probability, template count, and 
number of segments. This conclusion also follows from 
the similarity in form between Eqs. (3.13) and (3.311. 



To illustrate the relationship between p and p. Fig 
plots Zfa/^fa ^^'^ -^sV-^s as functions oi Pfd, and Fig, 
plots Pfd/pfd as a function of pfd and p/p. We see that 
both Zfa/zfa and N!./Ns are less than unity for pi^ > 0, 
while is generally greater than by a factor of ^ 4~ 
8. Note too that the right-hand side of Eq. (3.311 is 
divided by Rq w 0.56 [Eq. (|3.14|)]. This implies that p is 



dominated by contributions from linearly-polarised si, 
nals, as seen in the distribution of R?^^ plotted in Fig. U 
Therefore, we may think of p as estimating the sensitivity 
of a search to a population of constant-SNR signals (as 
for p), but where the signals are linearly polarised (hence 
the division by Rq), and where the search is performed 



with reduced false alarm threshold , 



■fa' 



a reduced number 
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FIG. 8. Ratios of the effective false alarm threshold Zf^/zfa 
(solid) and number of segments N!^/Ns (dashed) as functions 

of Pfd. 
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FIG. 10. Histogram of the relative difference between p 
given by Eq. (3.311, and the numerical solution to Eq. (3.9 1 



p, over the ranges 10"^^ < pt^^/Nt < 10~^ and 1 < iVs < 10 
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IV. ACCURACY OF ANALYTIC SENSITIVITY 
ESTIMATOR 



In this section, we validate the accuracy of the analytic 
sensitivity estimator p, derived in Section |II1D[ against 
the sensitivity calculated using two methods: numeri- 
cally solving Eq. (3.9), in Section IV A and performing 
software injections, in Section [lVB The validation is per- 
formed using a detection statistic with v — A, and for a 
target false dismissal probability of pfd = 0.1. 



A. Numerical solution of false dismissal equation 



FIG. 9. Contours of the effective false dismissal probability 
Pfd/Pfd ^ ^ function of pm and p/p. 



of segments iVg, and a greatly increased false dismissal 
probability p^^ . 



We note that Eqs. (3.25) and (3.28) impose restrictions 
on permissible values of pfd < 0.5 and p'^^ < \/e/2 « 0.82 
respectively. The latter restriction and Fig. [9] implies 
that, depending on the value of p/p, pf^ is further re- 
stricted to be less than 0.82/4 -0.82/ 8 w 0.1-0.2. In gen- 
eral, therefore, the use of Eq. (3.31) is restricted to val- 
ues of Pld 10-20%. In practice this is not an oner- 
ous restriction, as we are generally only interested in 
small false dismissal probabilities of either 10% or 5%. 
Gravitational-wave pulsar searches of LIGO and Virgo 
data (see Section |l] for references) have set upper limits 
with corresponding upper limit confidences of 90% and 
95% respectively. 



We first compare the sensitivity predicted by p against 
the sensitivity \/{f?) calculated by solving Eq. (3.9) nu- 
merically. We denote by p the value of \f{p^ which 
solves Eq. (3.9) for a target false dismissal probability 
Pid- We denote by ptd(/5triai) the result of computing 
Eq. (3.9) for a given trial value of p, denote d pt riai- For 
Ptriai = 0, Pfd(ptiiai) > Pfd, othcrwisc Eq. p?9l ) has no 
solution for the chosen pfd- We start by determining 
a Pmax such that Pfd(pmax) < Pfd, thus bracketing p 
to between and Pmax- We then use a simple bifur- 
cation search to converge to p, which terminates when 
the relative error between pjdifi) and the target pid is 
less than 10~^. Equation ( |3.9[ ) is solved for a grid of 
logarithmically-spaced values of pfa/iVt and N^- 30 val- 
ues in the range 10" < Pfa/A^t < 10"^, and 28 values 
in the range 1 < A^s < 10''. 

Figure [To] shows the relative difference between the an- 
alytic p, and the numerically calculated p. The maxi- 
mum relative error between p and p is 1.4% over the 
given range of Pfa/Nt and Ng. Note that p is consistently 
smaller than p, i.e. p slightly overestimates the sensitivity 
calculated using Eq. (3.9). 
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B. Software injections 

We next determine whether the sensitivity predicted 
by /5, and calculated numerically from Eq. (3.9) in the 



previous section, correctly predict the performance of a 
real gravitational-wave pulsar search pipeline. To do so, 
we perform software injection studies similar (with some 
simplications) to the Monte Carlo simulations used to set 
upper limits for gravitational-wave searches of LIGO and 
Virgo data (see Section |l] for references). We use software 
from the LA LSuit^ repository. 

First, the Makefakedata_v4 program is used to generate 
gravitational wave strain data, of timespan NgTs (see Ta- 
ble containing Gaussian noise (with a power spectral 
density of Sh = 1), and a simulated gravitational- wave 
pulsar signal, as it would be observed in the LIGO Liv- 
ingston detector. The strain amplitude ho is given by 
[c.f. Eq. 



ho 




(4.1) 



where Sh = ^,Ts is given values from Table [Tj and either p 
or p are calculated for pfd = 0.1 and values for pts^/Nt and 
Ns given in Table |l| Other parameters of the simulated 
signal are chosen uniformly from the following ranges: 
-1 < ? < 1, < V < 27r, < 00 < 27r, < a < 27r, 
-1 < sin (5 < 1, / = 100 Hz, and -10"* Hz s"^ < / < 
(with higher-order spindowns set to zero). 

Next, the ComputeFStatistic_v2 program, an implemen- 
tation of the J^-statistic, is used to perform a single- 
template search of the data generated by Makefake- 
data_v4, at precisely the sky position (a, S) and frequency 
evolution (/, /) of the simulated signal. (We consider the 
loss in sensitivity due to mismatch between the search 
template and signal in Section V C| ) A single value of 
the J^-statistic is returned by ComputeFStatistic_v2 and 
stored. For iVg > 1, since ComputeFStatistic_v2 cannot 
combine searches of multiple data segments, we instead 
run ComputeFStatistic_v2 iVg times on successive seg- 
ments, of timespans Tg, and add together the Ns returned 
J^-statistic values. Thus, the (summed) J^-statistic val- 
ues will follow distributions with Ngf = 4iVs degrees 
of freedom. 

Finally, we repeat the injection procedure 5000 times. 
We record the number of (summed) J^-statistic values 
which are below the false alarm threshold Sfa, as calcu- 
lated in Section [HI A| This fraction equals the false dis- 
missal probability as determined by the software injec- 
tions, which we denote Pfd,inj- If the sensitivity of Com- 
puteFStatistic_v2 is accurately estimated by p and/or p, 
Pfd.inj should be close to the target false dismissal prob- 
ability of pfd = 0.1. Different choices of pfa,/Nt are used 



Available from https : / /www. Isc-group .phys .UMin. edu/daswg/| 
P projects/lalsuite.htmll 



TABLE I. Validation of the analytic sensitivity estimator p, 
and the numerically-computed sensitivity p, using software 
injections. The injections are performed for three values of 
Pfa/Nt (listed in row 1), five combination of A^s and Ts (listed 
in columns 1 and 2), and once each using either p (top panel) 
or p (bottom panel). Values of Ap (or Ap) and Pfd.inj are 
given, for each of the three values of pta/Nt, in columns 3-4, 
5-6, and 7-8 respectively. 

N, Te Pfa/iVt = 10"' Pfa/iVt = IQ-' pta/Nt = 10'^° 

Analytic p 

Ap pfd.inj Ap Pfd,inj 

1 1 -5.% 13.2% -2.5% 12.4% -5.% 14.9% 
1 5 -2.4% 11.8% -2.8% 12.6% -2.9% 13.% 
1 10 -4.9% 13.2% -4.7% 14.% -2.5% 12.7% 
10 1 -1.3% 11.1% -1.7% 11.9% -1.7% 12.2% 
25 1 -0.16% 10.5% -0.93% 11.4% -1.5% 12.2% 
Numerical p 

Ap 

Pfd.inj 

Ap 

pfd.inj 

Ap 

Pfd.inj 



1 


1 


-3.% 


11.5% 


-4.1% 


12.9% 


-3.8% 


13.1% 


1 


5 


-4.7% 


12.4% 


-3.5% 


12.4% 


-3.4% 


12.8% 


1 


10 


-2.9% 


11.4% 


-3.% 


12.% 


-2.6% 


12.1% 


10 


1 


-0.2% 


10.1% 


-1.6% 


11.2% 


-1.8% 


11.5% 


25 


1 


0.36% 


9.8% 


-1.2% 


10.9% 


-0.99% 


10.9% 



to set different thresholds Sfa, and thus simulate the sen- 
sitivity of a wide-parameter-space search over A't tem- 
plates. 

Ideally, we would then refine the injected value of ho, 
increasing it if ptd,inj > Pfd, decreasing it if pfd,inj < Pfd, 
and then repeat the entire injection procedure, until 
Pfd,inj — Ptd- The rel ative error between the value of ho 



calculated using Eq. (4.1), and the value of ho arrived at 
by repeating the injection procedure, would then be equal 
to the relative error in the estimation of the sensitivity 
of ComputeFStatistic_v2. Since the injection procedure is 
time-consuming and computationally intensive, however, 
we instead re-compute p or p, as appropriate, using the 
false dismissal probability Pfd.inj determined by the in- 
jections. We denote by Ap = |p(pfd,inj) - p{P!d)\/p{Pid) 
the relative error between the value of p calculated using 
Pfd,inji and the value of p calculate using pf^; similarly 
for Ap. These quantities serve as a reasonable estimate 
of the error in the sensitivities estimated by p and p re- 
spectively. 

The results of the software injections are shown in Ta- 
ble [ij Both p and p underestimate, by < 5% and < 4.7% 
respectively, the ho required to achieve a false dismissal 
probability of 10%; except in one instance (last row of 
column 4), Pfd,inj > 10%. Nevertheless, these errors are 
still within the ^ 10% typical calibration error of gravi- 
tational wave detectors, e.g. LIGO [S7], and hence can be 
considered small. The difference between the predicted 
and actual sensitivity of ComputeFStatistic_v2 is likely 
because the J^-statistic values returned by ComputeFS- 
tatistic_v2 do not strictly follow a distribution; small 
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errors, on the order of a few percent, are introduced due 
to implementation details of the code [58l [59] ; see Sec- 
tion and Figs, [n] and [12] 



V. VALIDITY OF ASSUMPTIONS 

In the previous two sections, we have derived an ana- 
lytic estimator of the sensitivity of wide-parameter-space 
gravitational- wave pulsar searches (Section III), and val- 
idated its accuracy (Section |IVp . In doing so, we made 
certain assumptions: that the values of the detection 



0.2 



statistic returned by the search are statistically indepen- 
dent, that the network of gravitational wave detectors 
being searched is sensitive to gravitational waves from 
all sky locations, and that any signal present in the data 
is perfectly matched by at least one of the searched tem- 
plates. In this section, we investigate to what extent 
these assumptions are valid for real gravitational-wave 
pulsar search pipelines. 



A. Statistical independence of templates 

It was assumed, in deriving an expression for the false 
alarm probability in Section [ill A[ that the iVt values of 
the detection statistic s returned by the search are sta- 
tistically independent. This is not necessarily the case 
in practice. Gravitational-wave pulsar searches typically 
use template banks with small, e.g. 20% mismatches; 
templates nearby in parameter space will therefore have 
similarly-shaped waveforms. Matched filtering of nearby 
templates will therefore combine the same data with 
nearly-identical waveforms to produce corresponding val- 
ues of the detection statistic. There is the potential, 
therefore, for values of the detection statistic computed 
from nearby templates to be correlated with each other. 

To assess to what extent this effect is important, we 
perform four fully-coherent searches, using ComputeF- 
Statistic_v2 (see Section IV B I of 3 days of computer- 



generated Gaussian noise. The search parameter space 
is sky and frequency, with higher spindowns set to zero. 
Template banks were generated using the gridType=2 op- 
tion, which places templates over the sky using an adap- 
tive mesh. The first three searches are performed us- 
ing template banks with the following mismatches: an 
unrealistically-sniall mismatch of 1%, a realistic mis- 
match of 20%, and an unrealistically-high mismatch of 
500%. The fourth search repeats the first search using the 
1% mismatch template bank, with the following modifi- 
cation to the ComputeFStatistic_v2 program: before each 
value of the J^-statistic is computed, the input Gaussian 
noise data is regenerated, so that each J^-statistic is com- 
puted from independent instances of Gaussian noise. The 
frequency bands of the searches are chosen such that the 
searches return 4.2 x 10^ values of the J^-statistic. 

We partition the returned J^-statistic values into 2100 
blocks of = 2000 values, contiguous in frequency, and 
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FIG. 11. Histograms of the maximum values of the J^- 
statistic obtained in blocks of 2000 templates returned by 
an all-sky-frequency search of 3 days of computer-generated 
Gaussian noise. (Left) Histogram of the 1% mismatch tem- 
plate bank, computed from the same Gaussian noise (black), 
and from regenerated independent Gaussian noise (gray). The 
theoretical distribution is plotted for TV = 2000 (solid black 
curve), and A'^ = 400 (dashed black curve). (Right) His- 
tograms of the 20% mismatch template bank (black), and 
the 500% mismatch template bank (gray). The theoretical 
distribution is plotted for N = 2000 (solid black curve). 



select the maximum value of the J^-statistic in each block. 
If all the values of the J^-statistic are mutually indepen- 
dent, the distribution of the 2100 J^-statistic maxima is 
expected to be 33J 

p(2J-max) = Np{2T = 2J-„ax|2J- - x'[4]) 



p(2^<2J-„,,,|2^^x'[4]) 



(5.1) 



If, however, the values of the J^-statistic exhibit some mu- 
tual correlation, then we expect the distribution of the 



2100 J^-statistic maxima to be well-modeled by Eq. (5.1 1, 



but with an effective number of statistically independent 
templates N < 2000. This procedure was used by the 
gravitational-wave search for Cassiopeia A [12] to esti- 
mate the statistical correlation of the template bank. 

Figure [TT] plots histograms of the J^-statistic maxima 
obtained from the four searches. In the left-hand plot of 
Fig. 11 we see that the J-"-statistic values from the 1% 



mismatch search (without regenerating the input data) 
are highly correlated; they are best fitted by Eq. (5.1) 



with N — 400 effectively statistically-independent tem- 
plates. When the input data is regenerated before com- 
puting each J^-statistic value, the distribution is much 
closer to the expected distribution with N = 2000 statis- 
tically independent templates. This demonstrates that 
the origin of the correlations is that the same data are 
being used to compute the J^-statistic values of nearby 
templates. The small differences between the computed 
and expected distributions is likely due to the implemen- 
tation details of ComputeFStatistic_v2 (see Section IVB). 
In the right-hand plot of Fig. [TT] we see that the dis- 
tribution of the J^-statistic maxima, for both the 20% 
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FIG. 12. Histograms of the maximum values of the T- 
statistic obtained in blocks of 2000 templates returned by an 
frequency-spindown search of 7 days of computer-generated 
Gaussian noise. Details are the same as Fig. |11| 



and 500% mismatch template banks, are both close to 
the expected distribution with TV = 2000. This suggests 
that, while statistical correlation between J^-statistic val- 
ues is noticeable for very closely-spaced templates, it is 
less significant for realistic template bank mismatches. 

We repeat the above four searches, instead using 7 days 
of data and searching over frequency and spindown. Tem- 
plate banks were generated using the gridType=8 option 
to ComputeFStatistic_v2, which places templates over fre- 
quency and spindown using a lattice. The resulting dis- 
tributions of the J-"-statistic maxima, plotted in Fig. [12] 
are similar to those presented in Fig. [TT] except that 
the 20% mismatch distribution is less well-fitted by the 
N = 2000 distribution. The same template bank was 
used in the search for Cassiopeia A [T2j, which found 
the number of statistically-independent templates to be 
^ 88% of the total number of templates. 

While it is difficult to precisely quantify the effect of 
statistical correlations between templates on search sen- 
sitivity, for template banks with realistic mismatches the 
effect is likely to be small. As may be deduced from 
Fig. |6j a change in of e.g. 10% does not significantly 
alter the predicted p. Given that, as shown in Table |Ij 
both analytic and numerical sensitivity estimators accu- 
rately predict the sensitivity of a real gravitational-wave 
search pipeline to ~ 5%, it is reasonable to neglect the 
smaller effect due to statistical correlations between tem- 
plates. 



B. Isotropic sensitivity of detector network 

We now investigate whether it is reasonable to assume 
that a real gravitational-wave detector network, e.g. that 
of the LIGO and Virgo detectors, can be modeled by a 
network which is isotropically sensitive to gravitational 
waves arriving from all directions. While this assump- 
tion is reasonable for all-sky searches, since averaging 



over sky position is equivalent to averaging over detector 
orientation (see Appendix |A| , it may not be reasonable 
for searches targeting a single sky position. 

we numerically solved Eq. (3.91 for 



In Section IV A 



P — \/KP^)^ assuming an isotropically-sensitive detec- 
tor network, and averaging over the sky position (a, (5) 
and polarization angle i/j; we denote the values of p thus 
obtained by piso- We now solve Eq. (3.91, using the 
same algorithm detailed in Section |IV A for three de- 
tector networks: the LIGO Livingston detector, the two 
LIGO (Livingston and Hanford) detectors, and the three- 
detector network comprising the LIGO and Virgo detec- 
tors; detector locations and orientations are taken from 
|60j . We assume that all detectors in the network are 
equally sensitive, since we are concerned only with the 
effect of the geographic configuration of the network. We 
consider the following five search scenarios: 

1. A search using a long data segment length of 
Ts ^ days, covering the entire sky. For large Tg, 
the signal SNR becomes independent of a (see Ap- 
pendix [A|. We compute p averaged over the re- 
maining parameters, sin (5 and ■0. 

2. A search using a long data segment length of 
Tg S> days, and targeting a source at a known sky 
position but with an unknown polarization. We 
compute p for a grid of 41 linearly-spaced values of 
sin (5 in the range — 1 < sin 5 < 1, and average over 
^. 

3. A search using a long data segment length of 
Ts ^ days, and targeting a source with both a 
known sky position and polarization. We compute 
p for the same grid of sin 8 values as in scenario [2] 
and a grid of 10 linearly-spaced values of -0 in the 
range — 7r/4 < t/i < 7r/4. 

4. A search using a short data segment length of Tg = 
0.5 days, and targeting a source at a known sky 
position but with an unknown polarization. For Tg 
shorter than a day, the signal SNR is a function 
of a, sin 5, and ■0- We compute p for a grid of 20 
linearly-spaced values of a in the range < a < 
27r, 19 linearly-spaced values of sin (5 in the range 
— 1 < sin (5 < 1, and average over -0. 

5. A search using a short data segment length of Tg = 
0.5 days, and targeting a source with both a known 
sky position and polarization. We compute p for 
the same grids of a and sin 6 values as in scenario |4j 
and the same grid of values as in scenario [3j 

For each of the above search scenarios, we record the 
minimum, mean, standard deviation, and maximum of 
the relative error \p — pigo|/pigo and relative difference 
(p — pigo)/piso, over the grids of a, sin 5, and values 
given above, and over the grids of Pfa/A't and iVg values 
given in Section [IV A[ 
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TABLE II. Relative errors and differences in tfie sensitiv- 
ity estimated assuming an isotropically-sensitive detector net- 
work, under five search scenarios (see the text for details). For 
each detector network (column 1), the mean and standard de- 
viation of the relative error |p — piso|/piso (column 2), and the 
minimum and maximum difference (p — Piao)/piso (columns 
3 and 4) are listed. The abbreviations L, LH, and LHV 
indicate the LIGO Livingston detector, the LIGO detector 
network, and the LIGO- Virgo detector network respectively. 
Each block of the table corresponds to a search scenario, and 
a summary of each scenario is given just above each block. 

Network Error Ip - piso|/piso Difference (p - piso)/Piso 
mean ± stdv. minimum maximum 

1. Ts — >■ CX3; averaged sin 5, and 
L 0.0058±0.0037% -0.018% +0.017% 

LH 0.0062±0.0042% -0.02% +0.017% 
LHV 0.0063+0.0044% -0.02% +0.014% 





2. Ts 


— >■ oo; known sin 5; averaged tp 


L 


3.4 


+1.8% -4.5% +6.5% 


LH 


2.6 


+1.3% -4.9% +3.8% 


LHV 


3.7 


+2.% -7.5% +5.4% 




3. 


Ts — > cxj ; known sin 5 and 


L 


4.3 


+2.4% -11.% +6.5% 


LH 


2.8 


+2.% -4.9% +8.7% 


LHV 


3.8 


+2.4% -7.5% +9.3% 


4. 


T = 0.5 days; known a, sin 5; averaged tp 


L 


5.1 


+3.6% -6.5% +16.% 


LH 


4.7 


+3.% -10.% +10.% 


LHV 


4.1 


+2.3% -8.8% +7.4% 




5. T = 


0.5 days; known a, sin 5, and ip 


L 


6.3 


+4.7% -14.% +20.% 


LH 


5.6 


+3.7% -10.% +18.% 


LHV 


4.4 


+3.1% -9.9% +15.% 
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FIG. 13. Mean relative difference (p — piso)/piso in the sensi- 
tivity estimated assuming an isotropically-sensitive detector 
network, for search scenario [5] (see the text), as a function of 
sin (5, for the LIGO Livingston detector (gray), the LIGO de- 
tector network (black, dashed), and the LIGO- Virgo detector 
network (black). 
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FIG. 14. Mean relative difference (p — piso)/piso in the sensi- 
tivity estimated assuming an isotropically-sensitive detector 
network, for scenario [4] (see the text), as a function of a and 
sin S, for the LIGO- Virgo detector network. The sidereal time 
at Greenwich at the mid-point of the observation time is de- 
noted by (j)s- 



Tabic |TT] shows the mean and standard deviation of 
the relative error, and the minimum and maximum rel- 
ative difference obtained under the five search scenar- 
ios. For an all-sky search (scenario [T]) , the assumption 
of an isotropic detector network is an excellent one. For 
the remaining search scenarios, the error in assuming an 
isotropically-sensitive detector network increases as Tg is 
reduced, and as more parameters are set to fixed val- 
ues. While the mean error is limited to < 6.3% for all 
four scenarios, the maximum difference can be up to 20% 
for a single detector, although it reduces to 15% for a 
three-detector network. There are only a few (potential) 
gravitational-wave sources for which the polarization an- 
gle may be determined with any accuracy; noted exam- 
ples are the Crab and Vela pulsars, where X-ray obser- 
vations of the pulsar wind nebula provide information on 
the pulsar's orientation [101 Hi]. Thus, scenarios [sl an d [s] 
are less likely to arise in practice (but see SectionWl] for 
a discussion of these scenarios in relation to the Power- 
Flux upper limit procedure). We conclude that, while 
the error in assuming an isotropically-sensitive detector 
network may be acceptable in many cases, it can be sig- 
nificant for particular choices of fixed search parameters. 

13 plots the mean relative difference (p — 
^iso for scenario [2] as a function of sin 5, for the 
three detector networks considered. The difference be- 
tween p and piso is smallest at sin 5 ~ 0.5-0.55, the ap- 
proximate latitudes of the three detectors (~ 30-46°); a 



Figure 

Piso)/Pi 



signal originating from these declinations would therefore 
be located at the approximate maximum sensitivities of 
the detectors. For the same reason, a single detector at 
low latitude (LIGO Livingston) is more sensitive to sig- 
nals arriving at low declinations, explaining the increased 
sensitivity relative to that of an isotropic detector net- 
work. The addition of two detectors at higher latitudes 
(LIGO Hanford and Virgo) shifts the network configura- 
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FIG. 15. Histograms of the mismatch distribution of a 3- 
dimensional body-centered cubic lattice template bank (left), 
and an ad-hoc Gaussian mismatch distribution (right). 
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FIG. 16. Histograms of the relative difference between Pmism, 
calculated using the lattice template bank (left) and Gaussian 
mismatch distributions (right), and p (zero mismatch) over 
the ranges 10"^^ < Pfa/A^t < 10"^ and 1 < A^s < lO". 



tion towards improved sensitivity at higher declinations. 



Figure 14 plots the mean relative difference {p — 
Piso)/piso for scenario |4] as a function of a and sin 5, 
for the LIGO- Virgo network. Consistent with Fig. [TSj 
the network is more sensitive to signals arriving at the 
poles, and less sensitive to signals arriving at the equa- 
tor, than an isotropically-sensitive network. The change 
in sensitivity as a function of a and S is comparable to 
previous studies of gravitational-wave detector network 
configuration; see e.g. Fig 3 in jBT]. 



C. Template bank mismatch 

Finally, we consider the reduction in sensitivity due 
to mismatch between the searched templates and any 
signal which may be present. The mismatch fi is related 
to the difference between the SNR (p) of a signal with 
parameters p, and the SNR p^{p') recovered by a search 
template with mismatched parameters p' . It is given by 



p'{p)-p'{p') 
p^{p) 



(5.2) 



Substituting {p^)R for p{p) and (p^)-Rmism for pip'), 
we find that i?^;^^ = R^{\ — p); thus the eflfect of mis- 
match is to reduce by a factor \ — p. For a large tem- 
plate bank, p can be considered as a random variable 
drawn from a distribution p{p) characterizing the geo- 
metric arrangement of the templates in parameter space. 
Equation (3.9 1 is then modified to additionally average 



over mismatch: 



N^{p')R\l-p)]) 



a, sin 5,tp.^.fi 



(5.3) 



where (• • •)p = dppip) • • • , and Mmax is the maxi- 

mum allowed mismatch. 

We solve Eq. ( |5.3| numerically for praiem = \/Jp^, 
using a modified version of the algorithm described in 
Section |IV A| We compute /5,„ism for two examples of 
pip): the mismatch distribution for a template bank 
constructed using a 3-dimcnsional body-centered cubic 
lattice with /imax — 20%, from [33], and an ad hoc 
Gaussian mismatch distribution with a mean of 10% 
and standard deviation of 2%, restricted to the range 
< p < 20%. Histograms of these distributions are plot- 
ted in Fig. |15[ mismatch distributions for lattice tem- 
plate banks in other dimensions are plotted in |62| . The 
lattice template bank mismatch distribution has a mean 
mismatch of 10% and standard deviation of 4%. 

We find that taking mismatch into account increases 
Pmism? relative to the equivalent (zero-mismatch) p com- 
puted in Section IV A by on average 5.3 ± 0.10% for 
the lattice template bank mismatch distribution, and 
5.5 ± 0.10% for the Gaussian distribution. This ef- 
fect is of the same magnitude as the error in the an- 
alytic/numerical estimators compared to software in- 
jections (see Table |l]). We note that, by substituting 
Ro ^1 — (p) (where (p) is the m ean of the mismatch dis- 
tribution) for Rq in Eq. (3.31), we can predict the ob- 

reasonably accurately; 

This suggests 



served relative increase in p„ 
for (p) = 10%, p is increased by 5.4%. 
that it is reasonable to model the effect of template bank 
mismatch as a uniform reduction in recovered SNR (by 
■\/l — and that the mean of the mismatch distribu- 
tion is sufficient to quantify this reduction; this approach 



is also taken in [31|. Figure 16 plots histograms of the 
relative difference between pmism and p; the histograms 
are narrow, with no long tails, confirming that the most 
typical reduction in SNR is close to the average reduc- 
tion. 

Figure [4] plots the distribution of -Rfso.mism = -^?so(l ~ 
p), where p are drawn from the lattice template bank 
mismatch distribution. Relative to the zero-mismatch 



distribution, Rt 



iso,mism 



is reduced for circularly-polarized 
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signals (at high i?: 



-2 

isOjinism/ 



, and exhibits a broader peak 



at hnearly-polarized signals (at low R( 



iso,mism 



VI. DISCUSSION 

In this paper, we develop a new analytic method of es- 
timating the sensitivity of wide-parameter space searches 
for gravitational-wave pulsars. The new method avoids 
the unphysical assumption of constant-SNR signals, used 
by previous analytic sensitivity estimation methods, and 
hence can accurately (to < 5%) reproduce the sensitiv- 
ity estimated using Monte Carlo simulations, but with- 
out their computational cost. Additional effects, not 
included in the analytic model but which may be im- 
portant for real gravitational-wave pulsar searches, are 
investigated: the statistical correlation between values 
of the detection statistic due to closely-spaced template 
banks; the deviation of the directional sensitivity of 
real gravitational wave detector networks from an ideal, 
isotropically-sensitive network; and the loss in SNR due 
to mismatch between a signal and the searched tem- 
plates. While the first of these effects can be neglected, 
the remaining two may need to be accounted for where 
accurate sensitivity estimates are required. 

It is important to note that an accurate prediction of 
the absolute sensitivity of a search may not always be 
necessary. For studies of the optimal configuration of a 
hierarchical search (e.g. [271 EI]), the relative sensitiv- 
ity of different search schemes is usually more important. 
For these studies, use of the constant-SNR sensitivity 
estimator p may be acceptable; Fig. [T] shows that p re- 
produces to < 5% the correct scaling (i.e. that of p) of 
sensitivity with respect to pta/Nt and Ng. On the other 
hand, given the similarity between the expressions for 
p and p [compare Eqs. (3.13) and (3.31)], these studies 
may also be able to utilise the more accurate sensitivity 
scaling of p with little modification. 

This work has restricted its attention to detection 
statistics which follow distributions, which does not 
cover all gravitational-wave search methods. In partic- 
ular, the Hough semi-coherent method |26j computes 
a statistic, the number count, which is the number of 
coherently-analyzed segments where the coherent detec- 
tion statistic (which may either be power or the J-- 
statistic) exceeds a set threshold. It follows that the num- 
ber count follows a binomial distribution. Since, however, 
the Hough method is generally used to analyze a large 
number of data segments (i.e. Ns is large), the binomial 
distribution can generally be approximated by a normal 
distribution (see [HI), an approximation also used in Sec- 
tion |III| Thus, the analytic sensitivity estimation method 
developed here may be applicable to the Hough method 
with minor modifications. 

Searches performed using the PowerFlux method [HI 
[20l [22] employ a slightly different formula tion o f the fre- 
quentist procedure described in Section II C (To dis- 



ing discussion we refer to the upper limits produced by 
the frequentist procedure described in Section II C| as 
population-averaged upper limits.) An upper limit on 
ho is set, using the Feldman-Cousins method of confi- 
dence interval construction [63j . for every searched pa- 
rameter (typically sky position, frequency, and first-order 
spindown) , and assuming either linearly-polarized signals 
(^ — 0, with a range of tp) or a circularly-polarized signal 
(1^1 = 1). For each polarization case, the maximum value 
of ho over the searched parameters is chosen. Two upper 
limits are then quoted: a best-case upper limit, which 
assumes circular polarization; and a worst-cas^ upper 
limit, where ho is further maximized over the linear po- 
larizations. 

The PowerFlux upper limit procedure differs from the 
population-averaged upper limit procedure in two re- 
spects. (In the following comparison, we assume that the 
data being searched is relatively free of non-Gaussian in- 
strumental noise.) First, while the population-average 
procedure first selects the maximum value of the de- 
tection statistic s over the parameter space, and then 
computes an upper limit on /iq using s, the PowerFlux 
procedure first computes an upper limit on ho for ev- 
ery value of s, assuming a fixed polarisation, and then 
selects the maximum upper limit. By computing upper 
limits assuming a fixed sky position {a,S), corresponding 
to each value of s, and a fixed ip, PowerFlux is essentially 
performing search scenario [s] from Table |ll] (although 
with different search parameters, e.g. Tg — 30 minutes, 
Pid — 0.05). From the maximum relative difference in 
sensitivity quoted in Table [IT] for this scenario, we esti- 
mate that the effect of maximising over ho instead of over 
s raises PowerFlux upper limits by ~ 10-20% relative to 
population- averaged upper limits. 

Second, the population-average procedure produces a 
single upper limit averaged over all parameters; the Pow- 
erFlux procedure instead quotes upper limits for best- 
case (circular) and worst-case (linear) polarizations. By 
setting upper limits at fixed polarizations (i.e. fixed val- 
ues of ^), the PowerFlux procedure is setting upper limits 
on a population of signals with a fixed value of R^, i.e. a 
fixed SNR. We can therefore predict the upper limits on 
ho set by the PowerFlux procedure by replacing {p^)]^ 
with Pp.f.-RL(0 in Eq. it follows that Pp.f. is given 

by p/R\so{£,)i where p is given by Eq. (3.13). The best- 



and worst-case sensitivities are then given by 



PP.F. best 
PP.F. worst 



p 



R,U\i\ = 1) V 5 
-p 4 



(6.1) 
(6.2) 



Ratios of Power-Flux best- and worst-case sensitivities 
to the sensitivity predicted by p range over 0.39 < 



tinguish these two frequentist procedures, in the foUow- 



^ Actually, PowerFlux worst-case upper limits, as quoted in e.g. 
|14| . also account for the worst-case mismatch; in this discussion 
however we assume that no mismatch is present. 
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Pp.f. hest/p < 0.46, and 1.1 < Pp.f. worst//5 < 1.3, for 
the range of p/p plotted in Fig. [jj An injection study by 
Dergachev [M] found ratios consistent with these ranges: 
Pp.f. bost//5 ~ 0.43, and pp.p. worst/p 1.2. Further in- 
vestigation into the differences between the PowerFlux 
and population-averaged upper Umits methods would fa- 
cilitate direct comparison of the upper limits produced 
by different gravitational wave searches, e.g. between 

Eg and [ig. 

Recent searches for gravitational waves from known 
pulsars lini HI] use Bayesian inference to set upper lim- 
its. The frequentist method described in this paper con- 
structs confidence intervals which are derived from the 
probability of the data, given a particular hypothesis e.g. 
that a signal is or is not present in the data. In contrast, 
Bayesian inference directly calculates the probability of 
the hypothesis, given the particular data that were ob- 
served. Despite their very different interpretations, a 
study in [65 found that the frequentist and Bayesian 
methods produce similar upper limits, in the limit of 
large signal amplitudes. Further research is needed to un- 
derstand fully the relationship between these two meth- 
ods. 
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points toward zenith. The directions along which the in- 
terferometer is insensitive to any gravitational radiation 
are given by the unit vectors zLa{t) and ±6(t), where 



a{t) = 



ni{t) - n2(t) 
2sin(C/2) 



m = 



nijt) + n2{t) 
2cos(C/2) 



(A2) 



and C is the angle between the interferometer arms. Note 
that a{t) and b{t) are orthogonal, i.e. a{t) ■ b{t) = 0. As- 
suming that the gravitational wavelength is much larger 
than the interferometer arm-length, as is the case for 
ground-based detectors, the response matrix R is given 
in terms of these vectors by (e.g. [5B] ) 



a{t) (g) b{t) + b{t) (g) a{t)] 



(A3) 



where (8) is the vector outer product. The vectors can be 
expressed in terms of time-independent components: 

a{t) ~ 0,1 cos ils{t — to) + 0,2 sin fis(t — to) + ^3 i (A4) 

where the Earth rotates about the unit vector fig in a 
right-handed sense with angular frequency f2s, and ai = 
0.2 X Cls, fla = Os X ao, = (P's ■ and = a{t = 

to); similarly for b{t). 

The polarization matrices H+ and Hx may also be 
written in terms of certain vectors |58|,l67j. The directions 
along which a gravitational wave creates no space-time 
perturbation are given by the mutually orthogonal unit 
vectors ±£+ and ±y+, for a purely plus-polarized wave, 
and ±Xx and ±yx, for a purely cross-polarized wave. 
The two sets of vectors are related by 

*+\/2 = *x-yx, y+\/2 = *x+yx- (A5) 
The cross-polarisation vectors are given by 
Xx ~ ^cosijj + fjsmip , yx = "fj cos ip — ^ sin ip , (A6) 



Appendix A: Antenna-pattern functions and their 

averages 

We present expressions for the antenna-pattern func- 
tions F+ {t) and Fx (t) of an interferometric detector, and 
the averages of F+(t)'^ and Fx{t)^ over various parame- 
ters; see also the references in Section III A 



The antenna-pattern functions can be written in terms 
of a time-dependent detector response matrix R(t) and 
constant polarization matrices H_|_ and Hx: 

F+it)^tT{R{t)^H+), Fx(t) =tr(R(t)^Hx), (Al) 

where tr is the matrix trace and ^ denotes transposi- 
tion. We implicitly assume a coordinate system at rest 
with respect to the gravitational radiation, e.g. the Solar 
System bary center. 

Let ni{t) and n2{t) be unit vectors along the interfer- 
ometer's arms, such that the cross product ni{t) x fi2(t) 



where ^ — {n x f2g)/||n x ClsW, i] = ^ x n, and — n = 
:e_i_ X y+ = ^x X yx is the direction of propagation of 
the wave. The polarization matrices are then given by 



H = X g) y + y t 



) X . 



(A7) 



where we henceforth take F{t), H, x, y to mean either 
F+{t), H+,x+, y+ or Fx (t), Hx , Xx,yx , as appropriate. 



Combining Eqs. ( Al ), (TASl, and ( A7), we find that the 



antenna-pattern functions F can be written as: 



m 

sin^ 



[ait) ■ x] [b{t) ■ y] + [o{t) ■ y] [b{t) ■ x] . (A8) 



Note that when a{t), b{t) and x, y are parallel to each 
other (e.g. when a{t) is parallel to x, implying b{t) is 
parallel to y), F{t) achieves its maximal values of ±sinC. 

The averages of F{t)'^ over isotropic source sky posi- 
tion (a, S) and orientation ^, and over isotropic detector 
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location (<I>s, A) and orientation 7, are given by: 



where 



'a, sin S^ip 



= da d(sin(5) / 



2-K 



d(sinA) / d7 



dtp 



27T 



87r2 ' 



F^ 
8^ 



(A9) 
(AlO) 



where $s = ^s{t — to) is the local sidereal time at the de- 
tector. Note that these two equations can be transformed 
into each other if one makes the exchanges a O $s, 
(5 O A, and V' ^ 7- Note too that, since (i^^)Q,sin5,v av- 
erages over all possible orientations of x and y, it must be 
invariant to the orientation of a{t) and b{t); for the same 
reason, {F'^)^^,sinX.-y must be invariant to the orientation 
of X and y. Finally, note that F{t) is invariant if one 
exchanges a{t) and b{t) for x and y, and vice versa. We 
conclude that the averages (F^)a,sin5,V' ^iid (F^)*^, sin a, 7 
must be equal. To calculate e.g. {F"^) a, sin 5,ip, we choose 
a{t) — (1,0,0), b{t) = (0, 1,0) for convenience, and sub- 
stituting n — (cosacosJ, sinQ!COS(5, sin(5) into Eqs. (A6) 
obtain 



{-^ ) a, sill (5, ^/i — {-^ )<I>e!,sinA,7 



1 



■ sin^ C ■ 



(All) 



To compute the average of F{t)'^ over time, {F'^)t, we 
write d{t) = Ci(i)a'i, where ci{t) = cos51s(i — to), 

C2{t) — sinOs(< — to), C3{t) = 1, and to is chosen to be 
the mid-point of the observation time, at which ao is 
calculated; similarly for b{t). We can now write 



(F2),=sin2CE f dt[c.{t)c,{t)B,,]\ (A12) 



where we define 

Bij = (a,: • x) {bj ■ y) + (a^ • y) {bj ■ x) 



(A13) 



To write down the result of the integration in Eq. ( A12 1, 
we first define 



J±, = Br^' ± B,n^„ , (Al4a) 
K±, = B,,,„ ± B,„,, , (Al4b) 
S±^ = ^/B,,,„B,n,, ± B,,,,B,„,„ , (Al4c) 

where i = 1, 2, 3, and / denotes a permutation operation: 
1 = 3' = 2", 2 = 1' = 3", and 3 = 2' = 1". We also note 
the following identities: 



J_i — —Bii , J_2 — B: 



22 : 



^+3 — B, 



33 ■ 



Employing Eqs. (A14) and (A15I, wc find that 
(F2),=sin2C^C„sinc^, 

n=0 



(A15) 



(A16) 



Co = 
Ci = 



C2^'p-3J+3 + \{Kl,-Kl,), 
C3^l{J-3K+2-K+iK+3), 



Ca = ■^{Ji 



Ki 



(Al7a) 
(Al7b) 
(Al7c) 
(Al7d) 
(Al7e) 



For T ^ 1 sidereal day, (F^)^ sa Cq, which may also be 
written as 

Co = 9a|&|a;|yf[ + d,a\b\x\y\ 

+ 2(a|6i+ai&|-a|6j|) 

X {x\yl + xlvl - x\y\) , (A18) 

where 



a\\ = Jls • Qo , 

x\\ = rj, • £ , 



2a\ = l-a\. 
2x\ = \- x\ 



(A19) 
(A20) 



and similarly for 6||, y\\ and y^. Note that a||, 
Xjj, and y|| (and consequentially aj_, a;j_, and ?/_l) are 
independent of the Earth's angular displacement. We de- 
duce that Co cannot depend on parameters defined rela- 
tive to the Earth's angular displacement, which are the 
detector's local sidereal time <i?s, and the source's right 
ascension a. 



Appendix B: Limited inverse of the central ^ 
distribution 



The right-hand side of Eq. (3.4) is equal to the nor- 



malized upper incomplete gamma function: 

I I 2rAr r(7Vsl^/2,Sfa/2) 

We use a limiting form of the asymptotic inverse of the 
incomplete gamma function given in [5S], which is suffi- 
ciently accurate for small values of pfa/A't. It gives the 
following expression for Sfa: 



Sfa = NsiyX{T]) . 



where 



V = Vo + 



In 



X{x) = -W^i 



Nsvrjo \\{rjo)-l 
erfc-\2pfa/iVt) , 



exp —1 — 



(B2) 

(B3) 
(B4) 

, a: > , (B5) 
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and W-i is the —1 branch of the Lambert VF- function 
(see e.g [68]). For the purposes of this work, A(a;) is well 
approximated by the following two functions [see |69j . 
Eqs. (4.13.6) and (4.13.11)]: 



where g{x) = [1 — tanh5(x — 3)]/2. Equation (B2) is 
accurate to < 0.3% for values of pfa/Nt < 0.01. 



^^(^^ = ^ + ^+¥ + 36 - 270 
for X < 2y^, and 

\2ix)=y+{l + y-'+y-^)\ny, y=l + 



(B6) 



(B7) 



as a; — )■ cx). A sufhciently accurate piecewise expression 
for X{x) is 

Ai(x), x<2, 
g{x)Xi{x) + [l~g{x)]X2{x), 2<x<4, (B8) 
A2(x), x>4, 



As Ns — oo, 7/0 ^ and A(?7o) Xi{t]o) w 1 
With this approximation, rj w ryo, Sfa ~ A'si^(l + 



Vo) 



and the normalized false alarm threshold [Eq. (3.5)] is 
Zfs. ~ \/2erfc~^(2pta/A'^t). Thus, in the limit of large Ns, 
Zfa is independent of Ng. 
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